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Large spin systems can exhibit unconventional types of magnetic ordering different from the fer¬ 
romagnetic or Neel-like antiferromagnetic order commonly found in spin 1/2 systems. Spin-nematic 
phases, for instance, do not break time-reversal invariance and their magnetic order parameter is 
characterized by a second rank tensor with the symmetry of an ellipsoid. Here we show direct ex¬ 
perimental evidence for spin-nematic ordering in a spin-1 Bose-Einstein condensate of sodium atoms 
with antiferromagnetic interactions. In a mean field description this order is enforced by locking 
the relative phase between spin components. We reveal this mechanism by studying the spin noise 
after a spin rotation, which is shown to contain information hidden when looking only at averages. 
The method should be applicable to high spin systems in order to reveal complex magnetic phases. 

PACS numbers: 67.85.Fg,67.10.Fj 


I. INTRODUCTION 

Magnetic order in spin 1/2 systems is commonly asso¬ 
ciated with either a ferromagnetic phase or a Neel antifer- 
romagnet, depending on the sign of the exchange interac¬ 
tions. The situation is richer for spins greater than 1/2, 
and other types of magnetic order can arise at low tem¬ 
peratures. Spin 1 systems, for instance, can support spin 
nematic phases with vanishing average spin (i) [T]. The 
magnetic order is then characterized by a non-zero spin 
quadrupole tensor, Qij = \{siSj -\- sjSi) which deviates 
from isotropy even without applied held, i.e. it describes 
an object with the symmetries of an ellipsoid. In the sim¬ 
plest case, with axial symmetry, the spin quadrupole ten¬ 
sor has the same mathematical form as the orientational 
order parameter of nematic liquid crystals [5] . There is a 
preferred axis in space (the director) without a preferred 
direction along that axis. 

Spin nematic phases have been identihed in lattice spin 
1 models (see, e.g, [SHin]) or in spin 1 Bose-Einstein 
condensates (BECs) [H] with antiferromagnetic spin- 
exchange interactions [IMI]. In solid state systems, 
most magnetic probes couple only to the magnetization 
and are therefore unsuitable to reveal spin nematic or¬ 
der. In spin 1 condensates, equilibrium properties have 
been characterized by measuring the populations of each 
Zeeman state. This is not always sufficient to establish 
the nature of the magnetic order. For instance, in the 
so-called broken axisymmetry phase [52] , where all three 
Zeeman sublevels are populated, ferromagnetic or spin 
nematic behavior cannot be distinguished from the aver¬ 
age populations alone. 

In this article, we propose a method to reveal spin- 
nematic ordering (or possibly other types of unconven¬ 
tional magnetic order), and apply it experimentally to 
spin 1 atomic condensates. We show that the spin noise 
following a spin rotation contains information about the 
initial state, which can be retrieved with a suitable sta¬ 
tistical analysis. In spinor condensates, magnetic order 
follows from the emergence of a well-dehned phase rela- 
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FIG. 1: (Color online): (a): Sketch of the experimental setup. 
7500 Bose-condensed atoms are confined in a crossed 
optical dipole trap with a homogeneous static magnetic field 
along z. A resonant oscillating magnetic field along y drives a 
spin rotation of the initial equilibrium state, (b): Absorption 
image of the atomic cloud after Stern-Gerlach expansion in a 
magnetic field gradient, (c)-(e): Classical picture explain¬ 
ing the principle of our measurement, (c): The average spin 
(s) of the condensate created in a single realization can be 
decomposed into a longitudinal component = (sz) and a 
transverse component (s_l) = {sx)^7: + {sy)ey, the direction of 
which is given by the angle a. (d): From realization to real¬ 
ization, the angle a varies randomly while and |(si)| stay 
constant. The mean spin vector (s) thus samples a horizontal 
circle of radius |(sx)|. (e): A spin rotation of the initial state 
rotates this circle by an angle fit along the y axis. The fluc¬ 
tuations As/ after rotation are proportional to the squared 
radius of the circle through a simple geometrical relation. 


tion between the components of the spin wavefunction 
in the equilibrium state. This phase-locking mechanism 
is not caused by any external field, but emerges from 
the interactions between the atomic spins. We show ev¬ 
idence for such a mechanism in a condensate of spin 1 
^^Na atoms. 

The article is organized as follows. In Sectionjll] we 
recall results on the geometry of spin 1 wavefunctions. 









2 


which are used to give a quantitative definition of spin 
nematic order. We connect it to the standard treatment 
of spinor condensates at T = 0 , and discuss the effect 
of finite temperatures. In Section|lITj we describe the 
method used to extract informations about the magnetic 
order from a measurement of spin noise after a known 
spin rotation. In SectionjlVj we describe our experimen¬ 
tal apparatus and methods. SectionjV] describes a first 
analysis of our experimental results, where the fluctua¬ 
tions of magnetization after spin rotation are monitored. 
In Section|V^ another, more refined analysis is presented, 
where a maximum-likelihood estimation of the equilib¬ 
rium single particle density matrix is presented. Both 
methods reveal the underlying spin nematic character of 
the equilibrium state. Section [VlI| summarizes our find¬ 
ings. 

II. THEORETICAL DESCRIPTION OF 
ANTIFERROMAGNETIC SPINOR 
CONDENSATES 

The purposes of this Section are first, to give a precise 
definition of spin nematic phases in terms of spin observ¬ 
ables, and second, to connect this definition to experi¬ 
ments with spin I Bose-Einstein condensates at T = 0 
and at finite temperatures. We will assume here that 
the spin 1 bosons are confined in a state-independent 
trap, tight enough to prevent the formation of spin do¬ 
mains in the equilibrium state (single-mode approxima¬ 
tion) [23] . The condensate wavefunction is then given by 
the product of a spatial mode function (('(r), common to 
all Zeeman states, with a spin 1 wavefunction |^), which 
describes the internal degrees of freedom. An important 
feature of ultracold spinor gases is that the reduced (lon¬ 
gitudinal) magnetization, — n_i, is conserved 

by binary collisions driving the system to its equilibrium 
state [in HO]. Experimentally, we prepare a spin mixture 
well before the BEC forms in our evaporation sequence, 
allowing us to adjust the longitudinal magnetization 
between 0 and 1 (see Section[rv|) . 

A. Geometric description of spin 1 wavefnnctions 

We first give a more precise definition of spin nematic 
order, and connect this definition with spin observables. 
To that end, it is convenient to express a spin 1 state in 
terms of its components in the so-called Cartesian basis 
llx), I?/), |z)} [46j formed by the eigenstates of Sa with 
eigenvalue 0, where a = x,y, z. In this Section, we re¬ 
strict ourselves to the case of pure states for simplicity. 

A spin 1 state can be written in the Cartesian basis as 

[SlIlellSllEl] 

1^') = (ii-b w) • |r), (I) 

where the vectors u, v are real and obey u'^+v'^ = I. The 
vectors u and v are not uniquely defined. Performing a 


gauge transformation vk —> dt' = e*'>'tl' transforms u and 
V as u' = cos{'-f)u — sin( 7 )i; and v' = cos(j)v + sin( 7 )M. 
As a result, we can choose 7 such that u ■ v = 0 and 

ll«ll >11^11- 

The state of a spin 1 particle can be uniquely described 
by the average spin vector, (s) = 2uxv, and by the spin 
quadrupole tensor Qij = ^{siSj -|- SjSi) (TrQ = 2). In 
the cartesian basis, we have 

Qij — ^ij {'aiUj ViVj^^ ( 2 ) 

or in a more geometrical form, 

^ 1 — 1 -|- .A , , 

Q = — - — u(S>u-\ — v^v + w(S>W- (3) 

The orthogonal units vectors u = M/||it||, v = t)/||u||, 
w = (s)/ II (J) II define the eigenaxis of Q, with eigenvalues 
(1 — A)/2,(l + -A)/2 and 1. The alignment parameter 
A > 0, defined as A = 2 ||m|| — I, characterizes the 
anisotropy of spin fluctuations in the plane perpendicular 
to the mean spin vector. 

There are two simple limiting cases. The first one is the 
case of an aligned state (also called spin nematic or polar 
state in the context of spinor condensates mi where 
the spin wavefunction, Idt) = u • |r), is the eigenstate of 
s ■ u with eigenvalue zero. In such a state, the average 
spin vanishes, ||(s)|| = 0, and the spin quadrupole tensor 
is Q = 1 — u^u with eigenvalues 0,1,1. In the literature, 
it is common to call u the director field. The tensor Q, 
or equivalently the director u, plays the role of the order 
parameter for spin nematic states. 

The second limiting case is the one of an oriented 
or fully magnetized state, for which the average spin 
is maximal, |l(s)|| = I. This is achieved when ||m|| = 
||u|| = l/-\/ 2 j and also corresponds to a non-zero spin 

quadrupole tensor Q = ^ (1 — w) with eigenvalues 
1 / 2 , 1 / 2 , 1 . 

For a generic, partially magnetized state, one can 
quantify the proximity to one or the other limiting cases 
by the quantity A, which characterizes the amount of 
alignment present in a given state. For purely aligned 
states A = 1 while for purely oriented states A = Q. For 
a generic state, the alignment A and spin length |l(s)|| 
are related by 

{sf+A^ = l. (4) 

This shows that measuring the length of the mean spin 
vector (s)^ is fully equivalent to measuring the alignment 

A. 


B. Ground state of spinor condensates 

In the single mode approximation where atoms in dif¬ 
ferent spin states share the same spatial mode [55], we 
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FIG. 2: (Color online) Equilibrium population Uq (solid line), 
transverse spin length |l{sx)|l (dotted line) and alignment A 
(dashed line) of an antiferromagnetic spin 1 condensate versus 
longitudinal magnetization for a fixed value of q/Us = 0.2 

(solid lines). The critical magnetization separing the broken 
axisymmetry from the antiferromagnetic phase is ruz^c = 0 . 6 , 
marked by the vertical dashed line. 


parametrize the spin state of the condensate as 
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— np+rriz ^ife+a)/2^ 




l — np — rriz „i(&—a)/2 
2 ^ 


( 5 ) 


where 0 and a are relative phases [47j . The quantum 
state 1 C) corresponds to a mean spin vector (s) = 

(^x) (quantities in small letters are normalized by the 
total atom number TV). The mean transverse spin {s±) = 
{sx)^x + (sy)ey points in a direction determined by a and 
its length is determined by 0 , 


(sx)^ = 2no ^1 - no + \/(l “ - 'ml cos 0^ . (6) 

shows that the relative phase 0 also determines 
the alignment of the state |C). 

For a given magnetization ruz set by the preparation 
sequence, the equilibrium state |C) minimizes the spin 
mean field energy FImF) the sum of the spin-exchange 
interaction energy and of the quadratic Zeeman energy 
(QZE) energy in an applied magnetic field B [11], 



niz 


FIG. 3: (Color online) The phase diagram in the rriz — q 
plane, where the three sets of experiments we have performed 
are located, denoted as a,b,c (a: rriz = 0.33, q/h = 6.00 Hz; 
b\ ffiz — 0.73, q/h = 33.7 Hz; c: = 0.71, q/h = 3.84 

Hz;). In the gray area above the phase transition line both no 
and |{sx)| are nonzero, whereas both vanish below the phase 
transition in the zero temperature case. 


0 to 0 instead). This is equivalent to maximizing the 
alignment A introduced above. 

For a partially magnetized system with given mag¬ 
netization rUz, the competition between the two terms 
in Eq. ([^ drives a phase transition at a critical qc = 

Usd-ml) nKniiisiiiniiis]. At zero tempera¬ 
ture, the equilibrium population Uq is zero below qc (“an¬ 
tiferromagnetic phase”) and assumes a finite value above 
(“broken axisymmetry phase” ) muniiis], as illustrated 
in Fig.[^ Fig.shows the equilibrium population Uq, to¬ 
gether with the length ||(sx)|| of the transverse spin and 
the alignment A. Although the mean transverse spin is 
not zero above qc [see Eq. (|^], its value remains small 
because 0 stays locked to tt. As a result, the alignment 

A = no + a/( 1 - no)2 - ml, ( 8 ) 

which would reach 1 in the absence of other constraints 
(thus realizing pure spin nematic states), stays very close 
to the maximum value given the conservation of m^, 
A max — nA" — ml. This justifies using the transverse 
spin length to determine the amount of alignment present 
in the state |C), even when (s) d 0 - 


C. Finite temperatures 


Emf 
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A 
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(sx)^ - gno. 


( 7 ) 


up to terms that depend only on . For the experiments 
reported in this article, the interaction strength is Us/h « 
38 Hz (see Section [TV D[ ) and q/h « 4 Hz to 34 Hz. 

Antiferromagnetic interactions {Us > 0, the case of 
sodium atoms) favor minimizing the transverse spin 
length. According to Eq. §, this is achieved by lock¬ 
ing the relative phase 0 to tt independently of the value 
taken by no, a (ferromagnetic interactions would lock 


At finite temperatures, the T = 0 description of a 
spinor condensate should be modified in two ways. First, 
the spin state of the condensate is subject to thermal 
fluctuations, and second, the population of the conden¬ 
sate is thermally depleted. In this Section, we examine 
these two effects in order. 

We first discuss the thermal fluctuations of the spin 
state of the condensate, which is described by a finite 
temperature spin ensemble as studied in details in m- 
Close to the phase transition at gc, the population Uq > 0 
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FIG. 4: (Color online) (a) Partial condensed fractions 
for each Zeeman component mp = +1 (dashed red), mp = 0 
(dash doted green) and mp = — 1 (dotted blue line) and total 
condensed fraction (black solid). Here /N is 

normalized to the total number of atoms, (b) Transverse spin 
length (sj_)^ versus temperature. The calculation was done 
for a spherical trap of frequency u)/{2p) = 405 Hz, N = 7500 
atoms, fn^ = 0.33 and q/h = 6 Hz. In physical units, kpT = 
8to corresponds to T « 150 nK. 


which minimizes the free energy is small. The spin state 
of the condensate is then well described by a statistical 
mixture oi \N : 0 states, with an approximately Gaus¬ 
sian distribution of noi 0 m- 

We now discuss the thermal depletion of the conden¬ 
sate population. The single-mode approximation only 
describes the lowest energy “spatial mode” into which the 
atoms condense. Higher energy modes can be thermally 
populated, leading to a condensed fraction fc = Nc/N 
lower than one. Here N and W denote respectively 
the total number of atoms and of condensed atoms, irre¬ 
spective of their internal state. To describe the thermal 
component of the non-condensed cloud, we have adapted 
the Hartree-Fock (HF) description proposed in in 
the uniform case to our experimental situation (see Ap¬ 
pendix]^ for details). 

The results of this calculation are shown in Fig.|^ for 
parameters relevant to our experimental situation, where 
we plot the partial condensed fractions for each Zeeman 
component fc^^\ defined as the ratio of condensed atom 
number in state rrip to the total atom number. The 
condensed fraction in mp = 0 decreases first. Above 
kpT > I.SHlu, the mp = 0 component is purely normal 
and the condensate is formed hy mp = ±1 only. As found 
in |28j , the contribution of the thermal component to the 
average spin vector is oriented opposite to the average 


spin of the condensate. The total transverse spin is thus 
naturally reduced with increasing temperature |48| . In 
the regime we have investigated, the temperatures fulfill 
kpT q,Us- Asa result, the non-condensate spin vector 
is always much smaller in magnitude than its condensed 
counterpart, and we find that the main effect that reduces 
the length of the transverse spin vector is the reduction 
of the condensed fraction. The results of Section lllBI can 
be directly used, provided one replaces the total atom 
number N by the condensed atom number Nc < N and 
the reduced populations by their condensed coun¬ 
terparts. For a total condensed fraction /c = 0.8, (s_l)^ 
is reduced to about 57% of its zero temperature value. 


III. SPIN NOISE REVEALS SPIN-NEMATIC 
ORDER 

In contrast to the phase 0, which is locked to tt in 
equilibrium by the spin-exchange interactions, the phase 
a is expected to take random values from one realization 
to the next. When dealing with many realizations of 
the same experiment, the initial many-body state is thus 
characterized by a statistical mixture 

rather than a pure state \C^) with N bosons in the spin 
state 1C). Only three parameters {e.g., n-o, m^, (.s_l)^) 
are needed to characterize the ensemble, down from four 
to specify completely each member |C). In spite of the 
randomness of the spin orientation, these three param¬ 
eters can still be measured using spin rotation provided 
one goes beyond single-particle observables and measures 
spin noise (recent experiments used similar techniques to 
reveal squeezing [29H32]). 

Figure[J-e illustrates the method geometrically in 
terms of the mean spin vector (s). The mean spin vec¬ 
tor for a general spin 1 pure state |C) lies on or inside a 
sphere of radius one, with the phase a describing the az¬ 
imuthal angle of the transverse component of the mean 
spin vector (panel c). The ensemble of possible initial 
states with a uniform distribution for a lie on a circle of 
radius |(s_l)| around the z axis (panel d). In order to 
measure this radius, we rotate the state by a known an¬ 
gle nt around the y axis and measure the magnetization 
m), after rotation (panel e). As seen from the figure, the 
initial fluctuations of the transverse orientation map to 
fluctuations of m'^, which are readily measured. 

For a more quantitative description, we use the stan¬ 
dard angular momentum algebra to obtain the rotated 
operator 

s'^ = Rl{nt)sAy{nf) = cos{nt)s^ - sm{nt)s^. (lo) 

Here and in the following, primed variables denote quan¬ 
tities evaluated after the spin rotation is complete. We 
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now introduce a key assumption; the initial density ma¬ 
trix is invariant under rotation around the z axis. This 
is satisfied in particular by the density matrix in Eq. 
with a random phase a uniformly distributed in [ 0 , 47 r]. 
The value of an observable measured after averaging over 
many realizations of the experiment is 

= da{n^On), ( 11 ) 

where TZ = ^ symbol stands for 

a double average : the first one, denoted by ( • ), is the 
usual average over the quantum state before rotation for 
each realization, and the second one is done over random 
values of a arising from one experimental realization to 
the next. Defining an average in this way allows us to 
obtain formula expressing measurement results without 
specifying the initial state. 

Using this result, we find the average magnetization 
after the pulse, 

{m^)a = ^{S'z)a = cos{nt)mz, ( 12 ) 

which is independent of (s_l). However, the variance of 
the same quantity is given by 

Am'^ = cos^{m)Aml + siT?{m){Sl + S^), ( 13 ) 
where 

Si + sI=N + Nq + 2No (iV+i + 7V_i) 

+ -f h.c.) . (14) 

In other words, relying only on the randomness of a we 
find that the variance of the magnetization Am'^ after 
the pulse measures the initial transverse spin fluctua¬ 
tions. This result holds for a short enough pulse, such 
that one can neglect any other terms than the oscillating 
field in the Hamiltonian during the evolution time. 

It is convenient to rewrite the variance as 

^ sin^(Ht)(s_L)^ 

-I- cos^{nt)Aml + ^ sin^(Ht)As^, (15) 

with (sx)^ the squared length of the mean transverse 
spin, and with As^ = {Si + Sl)/N'^ — (sx)^ its vari¬ 
ance. For a spinor condensate with (sx) 7 ^ 0, the term 
on the first line dominates over the smaller noise terms, 
and Am'^ « ^ sin^(Ht)(sx)^. We thus expect that the 
variance Am'^ oscillates with the rotation angle and 
reaches its maximum for ^It = tt 12 where the slope of m(, 
versus is maximum. In our experiment, the last two 
noise terms in Eq. ( [I5| are typically dominated by the 
preparation noise on (which also introduces noise on 
no in the equilibrium state, and thus on (sx)). 


IV. EXPERIMENTAL TECHNIQUES 
A. Condensate preparation 

We prepare spinor condensates in a well-controlled ho¬ 
mogeneous static magnetic field B oriented along the 0 
axis [see Fig. We start from a precooled thermal 
cloud of ^^Na atoms in a crossed optical dipole trap [53] . 
The atomic cloud is partially magnetized, with a mag¬ 
netization friz ~ 0.5 on average resulting from previous 
cooling steps. We adjust the magnetization by either 
demagnetizing the atoms further with near-resonant RF- 
magnetic field sweeps, or by magnetizing it by evapora¬ 
tion in a magnetic field gradient (“spin distillation”) [2Q| . 
We are able to produce final magnetizations ranging from 
friz = 0 to friz = 1, with a typical error of 2 — 3 %. 

After preparing a spin mixture well above the critical 
temperature for Bose-Einstein condensation, the depth 
of the optical trap is lowered in a few seconds to perform 
evaporative cooling. A hold time of 3 s is added after 
the end of the ramp to ensure that the cloud reaches 
equilibrium [ 20 ] . At the end of the evaporation ramp, 
the atoms are confined in the crossings of the two beams 
of the dipole trap, where the trapping potential is well- 
approximated by a harmonic trap with average trap fre¬ 
quency a;/27r ^ 405 Hz (the trap frequencies are in the 
ratio 1 : 0.85 : 0.5). 

Experiments reported in this article are performed 
with “almost pure” Bose-Einstein condensates (BECs) 
containing typically 7500 atoms at a trap depth Vt/A:^ ~ 
400nK. By “almost pure”, we mean that no discernible 
thermal component can be observed in absorption im¬ 
ages. The measured condensed fraction /c = Nc/N is 
usually obtained by fitting a bimodal profile to absorp¬ 
tion images [34] . In our experiment, the contribution of 
the thermal component becomes difficult to detect for 
condensed fractions larger than « 0 . 8 , and the bi¬ 
modal fitting procedure unreliable. This sets a lower 
bound fc > 0.8 on the condensed fraction for the ex¬ 
periments presented in this article. 

We probe the sample using absorption imaging after 
free expansion in a magnetic field gradient, as shown in 
Fig.[^, and measure the normalized populations of 
each Zeeman component rriF = 0,±1 m- The three 
Zeeman components are imaged after releasing the cloud 
from the trap in the presence of a magnetic force sepa¬ 
rating the Zeeman components. Specifically, we apply a 
quadrupole field Bq = b'{2xex — yey — ze^) together with 
a uniform “separation” field with b' « 7G/cm 

and Bx « 3G. The resulting adiabatic magnetic po¬ 
tential is given by Umag = + .Bql « 

gF'mFfrB\Bx\ + gF'rriFyBb'x -!-•••, with gF = — 1/2 the 
Lande factor and with hb the Bohr magneton. The 
quadrupole and separation field are ramped up in a few 
milliseconds, while the bias field BBz applied during the 
experiment is simultaneously ramped down. 
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B. Experimental implementation of Rabi 
oscillations 


We apply a spin rotation using a radio-frequency (RF) 
magnetic field along y oscillating at the Larmor fre¬ 
quency. This RF field induces Rabi oscillations with 
Rabi frequency il. After a certain evolution time t which 
determines the rotation angle fit, we measure the final 
populations after spin rotation. The bias field is 
small enough to neglect the quadratic Zeeman shift {q < 
100 Hz) compared to the Rabi frequency (n/27r ^ 5 kHz). 
At the end of the pulse, the separation field B^ex is in¬ 
creased first, folllowed by the magnetic gradient used for 
SG imaging and by the decrease of the bias field 
The timing of the sequence is shown in Fig. Ramping 
up the separation field Bx is done with a linear ramp of 
T = 3 ms duration, sufficiently slowly to remain adia¬ 
batic with respect to spin flips {ujlT <C 1). The optical 
trap is switched off 10 ms after the end of the RF pulse 


(see Section IV C below). 


We have tested this sequence in two special cases, 
where all the atoms are initially in the mp = +1 state 
and or in the mp = 0 state. We are able to prepare these 
two states with little preparation noise, Aniz ^ 1%. The 
measured oscillations are presented in Fig.[^ The con¬ 
trast is close to 100 %, and we do not observe any sizeable 
dephasing of the oscillations after several Rabi periods. 
This shows that the assumption of adiabatic following 
when ramping up the different magnetic fields is valid. 


C. Influence of spin mixing after the spin rotation 

The sudden change of the spin state due to the spin 
rotation should in principle trigger a spin oscillation 
dynamics nzi I3SHSH] driven by spin-exchange interac¬ 
tions during the 10 ms hold time following the spin ro¬ 
tation. As seen before, the applied magnetic field is 
also changed after the spin rotation, from B = BzB^ to 
B — Bz^z Bx^x- The quadratic Zeeman energy q in¬ 
creases during this ramp, according to the curve shown in 
Fig.iD. This increase is fast compared to the time scale 
set by spin-exchange interactions, h/Us ~ 25 ms, and it 
reduces spin-mixing dynamics due to exchange collisions 
that would otherwise develop during the 10 ms hold time 
after the RF pulse. 

Nevertheless, a residual dynamics still takes place and 
modifies slightly the population no measured in SG imag¬ 
ing. Note that the effect of the spin interaction during 
the RF pulse is negligible {Us/hfl ~ 0.008). We model 
the spin-mixing oscillations using the theoretical frame¬ 
work given in [38] (see Appendix]^. An example for 
©i = TT is shown in Fig. The main changes in uq 
occur early in the ramp. Once q has settled at its final 
value qf ^ hx 2.5 kHz, the dynamics continue as a small 
amplitude oscillation of the population no around an off¬ 
set value (the so-called quadratic Zeeman regime [37)1. 
The oscillation amplitude is small (~ Ug/qf 0.015) 


and comparable to our detection noise. Ghanging the 
magnetic held to higher values would further reduce the 
amplitude without signihcantly changing the offset of Uq. 
Taking the long-time offset as the measured value of n-o, 
we hnd that the effect of the ramp amounts to increase 
the relative population in ng from its initial value by up 
to 0.05 for an initial angle 0^ = tt, a small but measurable 
change. 

We emphasize that the spin-mixing dynamics does not 
change the magnetization nriz of the system, but only the 
individual populations nmp- Therefore, the occurrence 
of spin mixing does not inhuence the analysis of the vari¬ 
ance of niz after spin rotation in SectionjV] On the other 
hand, it does affect the maximum likelihood analysis, as 
detailed further in section ED 


D. Determination of Us from spin-mixing 
dynamics 


We have measured directly the exchange interaction 
parameter Ug by deliberately inducing spin-mixing dy¬ 
namics and recording the oscillations of the normalized 
population uq after a sudden change (see Fig. Start¬ 
ing from a condensate with all atoms in the mp = 0 state, 
prepared as explained above at a bias field B « 282 mG 
[q/h Ki 22Hz], we first apply a spin rotation to produce 
a mixture with roughly balanced populations in all Zee- 
man states. This results in an initial state as given by 
Eq. 1^, with no « 0.38 and mz « 0. Spin-changing col¬ 
lisions produce high-contrast oscillations in the Zeeman 
populations, as observed in previous work for mz ^ 0 
[m na 131 [szi. The oscillation period has been pre¬ 
dicted analytically in [38], and is a function of no,mz,q, 
which are known, and of Ug, which is not. We extract 
Ug/h « 38 Hz from the measured period T^sc ~ 16 ms 
(see Fig.[^). 


V. SPIN NOISE MEASUREMENT OF 
SPIN-NEMATIC ORDER 


We now describe our experimental results on the mea¬ 
surement of the transverse spin using spin noise, as de¬ 
scribed in SectionUni In total we have taken three dif¬ 
ferent data sets for different initial magnetizations and 
magnetic fields which we label o, 6, c (see Fig.[^. The 
first two cases are above the T = 0 phase transition, 
while the third one is below. In each case, we drive Rabi 
oscillations with Rabi frequency O for an evolution time 
t, as described for quasi-pure spin states in Section |lVB[ 
and record the evolution of the relative populations 
after spin rotation. 
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t [ms] 


FIG. 5: (Color online) Rabi oscillation starting either from a state with all atoms in mp = +1 (a) or mp = 0 (b). The residual 
fluctuations are dominated by preparation noise, imperfections in the Rabi rotation parameters and the detection noise, all 
with roughly comparable contributions. 


A. Magnetization variance above the phase 
transition 


B. Spin thermometry 


We first focus on data set a. Fig.|^ shows typical raw 
data for the relative magnetization mi, (a) and the rela¬ 
tive population rig (b) for different rotation times t. As 
a result of the random orientation of the transverse spin 
(due to the random nature of a), large shot-to-shot fluc¬ 
tuations of the individual populations are observed. The 


mean magnetization behaves as predicted in Eq. (12). We 


extract the Rabi frequency from a cosine fit to the 
mean population m'^ (see Fig.|^). 

Fig.[^ shows the variance of mi,, displaying the ex¬ 
pected oscillations at twice the Larmor frequency. We 
compare the experimental results to the prediction of 
Eq. (blue solid line). The transverse spin length 


(sx)^ is computed with 0 = tt, with the measured friz 
and with the population rig found by minimizing Emf 
[49]. Eor comparison, we also show the transverse spin 
length for the same friz, but 0 = 0 (red dotted line) 
and for random 0 with uniform distribution (green dash- 
dotted line), that would correspond to a ferromagnetic 
system and to a non-interacting system (no phase lock¬ 
ing), respectively. Our measurements are best described 
by 0 = TT, as expected for antiferromagnetic systems 
in equilibrium. This shows that the system attempts to 
minimize its transverse spin, or equivalently maximize its 
alignment, thereby revealing spin nematic ordering. 

As seen from Fig.[^,c, data sets b, c show the same 
behavior as the case a discussed above, an oscillation of 
the variance with fixed amplitude. Data b is qualitatively 
comparable to a. Case c, taken below the T = 0 phase 
transition, deserves a separate discussion which we defer 
to Section lVCl 


We attribute the slight difference between the mea¬ 
sured amplitudes of the variance oscillations and the pre¬ 
diction of Eq. (15) for 0 = TT in Fig.|^ to a non-zero 
temperature. We addressed this point in details for data 
set a using the Hartree-Fock treatment of Section [Il C| 
Generally, we have found that increasing the tempera¬ 
ture reduces the transverse spin per atom. Experimen¬ 
tally, the condensed fraction can only be estimated as 
fc ^ 0-8 (see SectionIV). We show in Eig.[^ a shaded 
area where the lower limit corresponds to fc = 0.8 and 
the upper one to fc = 1, indicating that even a small 
non-condensed fraction leads to a measurable decrease of 
the oscillation amplitude. In fact, the oscillation vari¬ 
ance can be seen as a low-temperature thermometer. A 
temperature T « 80 nK (condensed fraction fc « 0.9) is 
found to reproduce the observed oscillation level (dashed 
line in Fig.[^). 


C. Magnetization variance below the phase 
transition 


For data set c, one would expect rio = 0 and (sx)^ = 0 
according to the T — 0 mean field picture. In contrast, 
we find a small initial population ug « 0.04, and an os¬ 
cillation of the magnetization variance with a small, but 
non-zero amplitude. The dotted lines in the figure corre¬ 
spond to the theoretical predictions which take the initial 
measured rig into account (corrected for the small shift 
in n g due to the spin changing collisions discussed in Sec- 
tion jlV C| and 0 = tt. 

A first explanation for this behavior could be the pres¬ 
ence of the thermal (uncondensed) component. In a 




























time after rotation pulse [ms] 


FIG. 6: (Color online) (a): Schematic diagram (not to scale) 
showing the experimental sequence. “RF” indicates the rf 
pulse inducing spin rotations, Bz is the bias held applied be¬ 
fore and during the spin rotation, and “SG” denote re¬ 
spectively the “separation held” and magnetic held gradient 
required for SG imaging, (b) The ramp of after the spin 
rotation results in a time-dependent Quadratic Zeeman en¬ 
ergy (QZE) q increasing within Sms after the end of the rf 
pulse (top panel). The evolution after the spin rotation of 
the normalized population no due to spin-mixing interactions 
calculated from Eq. (All is shown in the lower panel for an 
initial phase 0 ^ = tt, and an initial population chosen such 
that the hnal population is no « 0.43 (as measured for data 
set a), and Us/h = 38 Hz. 




FIG. 7: (Color online) Spin-mixing oscillations (a) and cal¬ 
culated oscillation period (b). A ht to a damped sinusoid 
is shown in (a) as solid line, and yields an oscillation period 
Tobc ~ 16 ms indicated by dashed lines in (b). 


tation of Am'^ = 0.005, comparable to the oscillation 
amplitude of the variance in Figure]^. 


VI. MAXIMUM LIKELIHOOD ESTIMATION 
OF THE DISTRIBUTION OF 0 


spinor BEG [39], spin excitations are phase-locked to 
the condensed components, and a finite transverse spin 
originating from the uncondensed component could con¬ 
tribute to our signal. However, from the Hartree-Fock 
calculations described in Section lHCl we found that the 
transverse spin of the uncondensed component remains 
very small for our typical parameters, and cannot explain 
the measured signal. 

A second explanation comes from a finite temperature 
of the initial spin state of the condensate, which is then 
described by a statistical ensemble rather than a pure 
state as described in m and Section [Trc| This leads to 
a finite population inmp = 0 even below the phase tran¬ 
sition. By numerically integrating the thermal distribu¬ 
tion described by the free energy given in [27] for a typ¬ 
ical temperature T = 80 nK, we find a finite population 
no = 0.016. This leads to a maximal variance after ro- 


A. Principle of the method 

We now turn to a more general statistical analysis 
based on maximum likelihood estimation (MLE), which 
allows us to estimate the distribution of the angle 0 in a 
more quantitative way. It takes all available data into ac¬ 
count, including the population Uq which was not used in 
the previous analysis. Given a set of measurements, the 
MLE method finds the most likely distribution among a 
set of parameter-dependent model distributions, thereby 
providing a statistical estimator for said parameters. 

We model the initial state by a density matrix 

= , 16 ) 

with an integration measure dQ = driqdmzdQda. We 
assume for simplicity that the probability density func- 
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FIG. 8: (Color online) (a): Magnetization mi, and (b): relative population n'o in the mp = 0 state versus duration of the 
Rabi pulse (or equivalently, rotation angle). The small blue dots correspond to single-shot measurements, while the larger red 
circles correspond to the average (mi,,nQ) over all measurements for each pulse duration. The solid line shows a cosine fit to 
the average mi,, from which we extract the Rabi frequency fl. (c): Variance of mi, (blue circles) versus duration of the Rabi 
pulse oscillating at t wice the Rabi frequency. The blue solid line corresponds to the theoretical prediction at zero temperature 
obtained from Eqs. ( 6|15 1 and an initial phase 0 = tt. The red dashed line (0 = 0) and green dash-dotted line (random 0) are 
shown for illustrative purposes. The data are from set a, with m^ = 0.33 and q/h = 6.0 Hz {B = 147mG). 



t [ms[ 


FIG. 9: (Color online) Close-up view of the magnetization 
variance for data sets a, b and c. In a and b, the solid blue line 
is the zero temperature theory for an initial angle of 0 = tt 
( antiferromagnetic interactions). In a, the shaded area corre¬ 
sponds to the prediction of our Hartree-Fock model at finite 
temperature assuming a condensed fraction of fc > 80%. The 
dashed line gives the best agreement corresponding to a tem¬ 
perature T ~ 80 nK. In c, the data are shown for parameters 
below the T = 0 phase transition line. The dotted curve 
is the theoretical expectation from Eq. (151 taking the initial 
population no into account. 


tions Ginojiriz) and P{Q) are Gaussians. We note that 
the equilibrium density matrix of a finite-temperature 
spin ensemble is well-approximated by Eq. (16) with a 


Gaussian weight function m- The joint probability 
density G{nQ,mz) is peaked around the average value 
(nQjfriz) with Uq the population minimizing Emf, with a 
finite width mostly due to experimental imperfections in 
the preparation sequence. The covariance matrix char¬ 
acterizing G{nQ,mz) is extracted from the experimen¬ 
tal data. At r = 0, P{Q) is given by a Dirac delta, 
P(0) oc (5(0 — tt) but acquires a finite width at finite T 
(see SectionVID below). The mean value 0 and stan¬ 


dard deviation ere of P{Q) are the unknown parameters 
to be estimated. Due to the periodic nature of 0, our 
choice is sensible only when P( 0 ) is peaked around the 
mean, i.e. gq ^ 27r. 

We use a Monte Garlo method to sample the initial 
distribution in Eq. (16|. Eor a given 0 and a mea¬ 


surement time ti (rotation angle Glti), the initial state 
(no,TO 2 ,Q(, 0 ) is propagated in time using the rotation 
operator. Here we assume that the spin rotation is per¬ 
fectly known, with rotation axis y and a rotation angle 
extracted from the fit to as before. Spin-mixing 

dynamics just after the spin rotation slightly change 
the relative population tiq, and is taken into account 
in the propagation. After convolution of the final re¬ 
sults with our known measurement noise, we get a condi¬ 
tional probability density p*. (m(,, ng| 0 , cto) for the mea¬ 
sured (nQ,m(,). Given a set of independent observations 
{m(, j,ng j}, we can construct a (log) likelihood function 

logL(0,cre) = ^logpt,(mW,noJ0,cre)- (17) 


The distribution that accounts best for the observed re- 
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suits is found by maximizing this function. 

Since the estimator strongly depends on the chosen 
probabilistic model, it is important for this model to be 
close to the physical reality. In the following we motivate 
the model used in the MLE before discussing the results. 


B. Model for the initial distribution 



FIG. 10: (Color online) Initial distribution C/(no,mz) for data 
set a. (a) Measured initial populations are indicated by black 
circles. The color shows the two-dimensional histogram of the 
simulated 10® initial points used in the Monte Carlo method. 
These points are drawn from a two-dimensional Gaussian dis¬ 
tribution estimated from the initial measurements, and sam¬ 
pled on a square grid with a step size 2.5 x 10“®. The dashed 
lines limit the allowed area (|mz| < 1 — no), (b): Marginal 
histograms of measured initial populations for m^, and (c): 
no. The red line indicates the marginal distributions used for 
the Monte Carlo analysis. 


The distribution of initial states is probabilistic due to 
three different effects. The first effect is intrinsic to our 
theoretical model where the initial angle a takes ran¬ 
dom values from one realization to the next. The second 
probabilistic effect is due to experimental imperfections, 
mainly fluctuations of (from the preparation process 
and the subsequent evaporation), or fluctuations in the 
spin-spin interaction energy 14 (due to fluctuations of the 
total atom number or of the confinement strength). Such 
fluctuations result in correlated fluctuations in no due to 
the system exploring different minima of the mean held 
energy. We stress that the marginal distribution P(0) 
is a priori not affected by these huctuations. A third 
random element originates from the hnite spin tempera¬ 
ture as described in Section [lI C| which allows the system 
to explore states situated away from the minimum. The 
second and third effect are more pronounced close to the 
phase transition m- 

We hnd empirically that the initial joint distribution 


of no and rn^ in Eq. (16) is well described by a two- 
dimensional Gaussian G{no, rriz). The mean and covari¬ 
ance matrix characterizing G are calculated from the 
measured data without spin rotation. We account for the 
spin changing collisions discussed in Section [lV C[ which 
affect the measured “initial distribution”, i.e. the distri¬ 
bution observed without any spin rotation. Specifically, 
for each values of 0, and no, the mean field equations 
(A1) are used to find the initial value noy that leads to 
the measured one, no{t = 0). The known values of q and 
the measured value of Ug are used as fixed inputs for this 
calculation. The initial distribution t/(no,i, Wz) deduced 
in this way is shown in Fig.[^ We estimate that exper¬ 
imental imperfections dominate the initial distribution 
G{no,m^). 


C. Monte Carlo approach 


To compute the evolution of a given initial state un¬ 
der spin rotation, we use a Monte Carlo approach. The 
initial density operator is sampled by drawing random 
numbers (no,mz,Q;) according to our assumed probabil¬ 
ity distributions (see Figure [l0| and assuming a certain 
value for 0. This determines an initial mean field state 
\C^)- Using the known evolution under spin rotations, 
we propagate this state in time for a given ti to arrive 
at the final mean outcome populations (nQ,m(,) as the 
expectation values of the corresponding operators in the 
time-evolved mean field state. In our numerical imple¬ 
mentation we use a typical number of ^ 10® Monte Carlo 
samples to reconstruct the final statistical distribution of 
the measurement outcomes. Spin mixing collisions as 
discussed in Section lTV Cl are also taken into account to 
obtain the final simulated distributions. In the Monte- 
Carlo simulation, the spin state found after rotation is 
used as initial condition to solve the mean field equa¬ 
tions (Al) describing the spin dynamics. We arrive in 
this way at a distribution of Uq corrected for the effect of 
spin changing collisions, typically by a few percents. 


We evaluate the final populations for each realization 
using expectation values. Doing so, we neglect the effect 
of quantum fluctuations on the final results, which are on 
the order 1 / and small for our typical atom num¬ 

bers of particles {Nmp ~ a few thousands) when com¬ 
pared to the noise level of our population measurements. 
The measurement noise, caused by a combination of pho¬ 
ton shot noise and small spatial intensity fluctuations of 
the laser pulse used for absorption imaging, is typically 
Arimp ~ 1% for the normalized population in Zeeman 
state nip. We include this noise in our model by convolv¬ 
ing the simulated measurement outcome by a Gaussian 
distribution. This leads to a conditional probability den¬ 
sity pt;(nQ,m(,|0) for the measurement outcome which 
depends on the initial phase 0, which is then multiplied 
by the distribution P(0) to obtain pt-{nQ,m'^\Q,aQ). 
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FIG. 11: (Color online) Comparison for data set a of the Monte Carlo simulated populations with measured data (black circles) 
after Rabi-rotation for different assumed initial angles 0. The three panels show the color-coded two-dimensional histograms 
of the Monte Carlo simulations for 0 = 0 (a), 7 r/2 (b) and n (c). Best agreement, i.e. the maximum likelihood is found for 
0 = TT. The example is taken for Qt ~ 7r/2 where the sensitivity is the highest. 


D. Results of the MLE 


VII. CONCLUSION 


We model the distribution P(0) by a truncated Gaus¬ 
sian with a mean value 0 and a standard deviation 
tJe- For the three data sets a, b and c, the log like¬ 
lihood is shown in Fig.[l^ versus (0,(Te)- The maxima, 
shown in Fig.[^ are found for (0, cto) = (I.OItt, O.OSStt), 
(O.SOtt, 0 . 3477 ryand (I.OStt, 0.2107r) for data sets a, b and 
c, respectively. These results are in full agreement with 
the conclusion drawn from the variance analysis, confirm¬ 
ing the locking to tt of the relative phase 0 . 


In all instances, the MLE is maximum for a finite width 
P(0) which is not expected in the standard T = 0 de¬ 
scription. We conclude this Section by discussing possible 
explanations. First, it may be caused by an underestima¬ 
tion of the noise sources in the system. As seen before, 
the probability distribution pt-{nQ,m'JQ) is almost sym¬ 
metric in 0 with respect to 0 and tt. The presence of 
fluctuations (induced for example by experimental im¬ 
perfections) not included in our model always bias the 
estimator away from 0 = tt. We thus infer that underes¬ 
timated or unconsidered noise in our probabilistic model 
will result in a broadening of the estimated distribution 
P(0). A second, more fundamental effect comes from the 
finite temperature of the initial spin ensemble (see Sec¬ 
tioning. The marginal distribution of 0 obtained nu¬ 
merically m is a bell-shape curve centered at tt, reason¬ 
ably approximated by a Gaussian with root-mean-square 
(rms) width « ^ksT/NUs- Using T « 80nK and the 
experimental parameters of data set a, we find a width 
« 0.1 comparable to the results of the MLE. 


In conclusion, we have shown the existence of spin- 
nematic ordering in antiferromagnetic spin 1 BECs, or 
equivalently of a phase locking between the Zeeman com¬ 
ponents caused by spin-exchange interactions in the equi¬ 
librium state. Our experimental method combines spin 
rotations with a statistical analysis, either based on the 
spin moments or on a maximum-likelihood estimation of 
the probability density function characterizing the ini¬ 
tial spin state of the condensate. Our method is not 
restricted to single-mode condensates or to spin 1 atoms, 
and could be used to reveal other types of spin order¬ 
ing. We remark in particular that measuring the spin 
variance gives access to a quantity (the squared trans¬ 
verse spin length) which can be used to characterize other 
phases than a fully condensed state. The expression of 
the transverse spin operator in Eq. (14) shows that mea¬ 
suring the spin variance gives access to the “spin singlet 


amplitude” [iQlIlT], which appears in stud¬ 

ies of fluctuating systems beyond mean field (spin liquid 
in one dimension |42j, or spin-singlet Mott states in op¬ 
tical lattices, for instance El)- 


Acknowledgments 


We acknowledge support from IFRAF, from DARPA 
(OLE program), from the Hamburg Center for Ultrafast 
Imaging and from the ERG (Synergy grant UQUAM). 












12 



O/tt O/tt 0/7r 


FIG. 12: (Color online) Results of the maximum likelihood determination of P{Q), expressed in terms of the mean value 0 
and standard deviation uq for the datasets a, b and c. In all cases the maximum, i.e. the estimated phase 0, is found close to 
the theoretical predicted value of tt (black cross). The contour lines delineate the 67%, 95% and 99% confidence area. 
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FIG. 13: (Color online) Most likely values of 0 found by the 
MLE algorithm for data set a, b, c. Error bars indicate the 
67% confidence bounds. 


Appendix A: Calculation of spin-mixing dynamics 


To quantify the impact of spin-mixing oscillations on 
the measured no, we use the theoretical framework given 
in [3H] . The evolution of an initial state of the form given 
in is described by the two Josephson-like equa¬ 

tions , 


= 2UsnoyJ (1 - noY - ml sin(0), (Al) 


= -2q(t) + 2Us{l - 2no) 


(A2) 


, ^rr (1 - ’^o)(l - 2no) - ml 

= +2C^s'-=—cos(0), 

\/(l - noY - ml 


with no(0) = no,i and 0(0) = 0^. We solve Eqs. (Al) 
numerically with q{t) as shown in Fig.[^ to compute the 
evolution of no- 


Appendix B: Hartree-Fock model of a spin 1 gas at 
finite-temperatures 

The model of [2H] treats the non-condensed cloud 
as a gas of non-interacting free particles evolving in a 
self-consistent mean field potential accounting for spin- 
exchange interactions |28j . Importantly, this mean field 
potential is not diagonal in the Zeeman basis due to 
spin-mixing interactions. The thermal component can 
in principle develop non-zero coherences due to interac¬ 
tions with the condensate and therefore a non-zero aver¬ 
age spin. The quantity of interest is the single-particle 
density matrix, 

(Bl) 

with 4> the condensate wavefunction, with pm}n the con¬ 
tribution of the thermal component, and where m,n = 
0, ±1. The density in each Zeeman component m is de¬ 
termined by the diagonal terms pmlm and the transverse 
spin by the off-diagonal coherences 

With respect to the full HE model laid out in [^, we 
make two additional simplifying assumptions. First, we 
assume that the single-mode approximation holds for the 
condensate wavefunction [^. This amounts to setting 
4>m{r) = •\/Ac0(r)^m, as done in the main text. The 
single mode wavefunction (/) determining the condensate 
spatial distribution is computed numerically by solving 
the GP equation 

2MNa 

with MNa the mass of Sodium atoms. The spinor part 
C,m is found from the single-mode theory using Ug = 
NcPs J r\(j}\^. The coupling constants g,gs are pro¬ 
portional to the scattering lengths a « 2.79 nm and 
as « 0.1 nm |33] with a proportionality factor 47rfi,^/MNa- 
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Second, we neglect the contribution of the thermal cloud 
to the mean-field potential (“semi-ideal” model |44j). Far 
from Tc, this is expected to be an accurate approximation 
|45j . Finally, we perform the calculations for a spherical 
trap. Although the trapping potential used in the exper¬ 
iment is not exactly isotropic, we do not expect that this 
affects strongly the results (in the Thomas-Fermi regime, 
for instance, only the average trap frequency matters to 
compute thermodynamic quantities |45j ). 

The excitations modes and energies E^, are solu¬ 
tions of the eigenproblem 

. +V{r) + A(r)) (B3) 


where the matrix A, explicitely given in |28j . depends 
on the condensate wavefunction (p{r) and on g,gs- Di¬ 
agonalizing this equation, we obtain the single-particle 
density matrix of the thermal component as 


ul^Hr)NBE{E^) (B4) 


with Nbe{E) = — 1) the occupation number 

for each mode v. 
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